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ABSTRACT 


2  Two-dimensional  idealized  dry  and  moist  numerical  simulations  are  performed  and  analyzed 
with  a  nonhydrostatic,  fully  compressible  spectral  element  model.  The  dry  numerical  tests 
4  consist  of  a  linear  hydrostatic  mountain  wave  and  a  squall  line  is  the  basis  for  the  moist 
simulations.  A  desired  spatial  accuracy  with  a  spectral  element  model  is  determined  by  two 
e  parameters,  a  number  of  elements  ( h )  and  a  polynomial  order  ( p )  of  the  basis  functions.  In 
this  paper,  the  range  of  average  horizontal  resolution  varies  from  0.2  km  to  10  km. 
s  Dry  experiments  are  compared  to  an  analytic  solution  for  accuracy.  It  is  found  that  the 
nominal  resolution  ( Ax )  of  less  than  2  km  is  sufficient  to  minimize  the  error,  while  resolutions 
io  of  500  m  or  less  show  no  additional  error  reduction  and  are  computationally  expensive.  When 
compared  at  a  similar  spatial  resolution,  the  computational  cost  of  the  spectral  element  model 
12  compared  to  a  finite-difference  model  is  an  order  of  magnitude  larger,  but  the  accuracy  gain 
is  significant  with  the  error  an  order  of  magnitude  smaller  for  the  spectral  element  model 
14  when  Ax  is  less  than  1  km.  If  the  acceptable  error  is  known  a  priori  the  spectral  element 
model  is  less  costly  compared  to  the  finite-difference  model. 
i6  Evolution  of  a  simulated  squall  line  is  compared  across  the  h  p  space  and  evaluated 
with  the  help  of  three  integrated  quantities:  total  precipitation  accumulation,  maximum 
is  vertical  velocity  and  maximum  precipitaton  rate.  The  squall  line  is  adequately  resolved 
when  the  nominal  resolution  (Arc)  is  less  than  2  km,  but  in  addition,  the  polynomial  order 
20  ( p )  needs  to  be  at  least  5.  The  analysis  of  the  integrated  quantities  across  the  parameter 
space  consistently  shows  a  gradient  with  respect  to  h  at  a  fixed  p  value  (e.g.  less  rainfall, 
22  stronger  maximum  vertical  velocities  and  weaker  maximum  rain  rates  with  increasing  /?,). 

The  nonlinear  nature  of  moist  processes  is  responsible  for  this  resolution  dependence  as  a 
24  result  of  localized  buoyancy  sources,  evident  in  spectral  analysis  of  the  time-  and  height- 
averaged  vertical  velocity. 
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1.  Introduction 


Numerical  models  used  for  mesoscale  weather  forecasting  can  be  assembled  into  two 
groups  depending  on  the  approach  used  to  solve  the  governing  Navier-Stokes  equations. 
In  the  first  group,  the  equations  are  kept  in  the  differential  form  and  both  temporal  and 
spatial  derivatives  are  approximated  using  finite-differences  (e.g.  COAMPS,  Hodur  (1997), 
WRF,  Skamarock  et  al.  (2005),  MM5,  Dudhia  (1993),  MC2,  Benoit  et  al.  (1997),  LM,  Doms 
and  Schattler  (1997)).  In  the  second  group  are  methods  based  on  an  integral  form,  less 
frequently  used  in  mesoscale  forecasting,  which  includes  spectral  (Aladin,  Bubnova  et  al. 
(1995)),  pseudospectral,  finite-element,  spectral  element  and  finite-volume  models. 

Spectral  element  models  have  been  traditionally  used  in  computational  fluid  dynamics 
and  more  recently  also  in  computational  geophysical  fluid  dynamics.  For  example,  atmo¬ 
spheric  phenomena  have  been  studied  on  the  global  scale  (e.g.  Giraldo  and  Rosmond  (2004); 
Fournier  et  al.  (2004)),  and  more  recently  also  on  the  mesoscale  (Giraldo  and  Restelli  2008; 
Giraldo  et  al.  2010).  Examples  for  the  ocean  include  a  density  current  model  on  a  local  scale 
(Ozgokmen  et  al.  2004),  and  general  circulation  model  (Dupont  and  Lin  2004;  Curchitser 
et  al.  1998). 

Advances  in  high  performance  computing  have  led  to  substantial  increase  in  the  number 
of  computational  cores  and  to  a  smaller  extent  improved  speed  per  core.  Availability  of 
computational  resources  is  facilitating  horizontal  resolution  refinement  for  the  numerical 
weather  prediction  models  on  the  global  scale.  In  the  future,  one  can  foresee  a  natural 
merging  of  global  and  local  area  modelling  efforts.  Ideally,  future  unified  models  will  have 
a  flexibility  to  allow  for  varying  resolution  with  grid  refinements  over  areas  of  interest.  In 
addition,  the  model  should  be  able  to  utilize  many  cores  and  also  scale  well. 

The  spectral  element  (SE)  method  has  the  potential  to  meet  this  changing  and  chal¬ 
lenging  computational  paradigm.  This  method  provides  a  flexible  platform,  which  supports 
unstructured  grids  and  provides  flexibility  to  adjust  the  accuracy  of  the  dynamical  core  with 
a  simple  change  of  a  control  parameter.  Moreover,  the  communication  requirements  in  par- 


allel  processing  are  minimized  because  the  adjacent  elements  share  only  the  boundary  points 
2  with  no  additional  interior  points  to  be  exchanged  (Kelly  and  Giraldo  2011).  An  additional 
novel  characteristic  of  the  SE  model  presented  in  this  paper  is  the  vertical  discretization, 
4  which  is  traditionally  based  on  either  a  finite-element  (Beland  and  Beaudoin  1985)  or  finite- 
difference  formulation  (Kim  et  al.  2008).  The  model  applied  in  this  study  is  two-dimensional 
e  and  spectral  element  in  both  horizontal  and  vertical  direction.  To  our  knowledge,  this  is  the 
first  fully  compressible,  nonhydrostatic  spectral  element  model  which  also  includes  cloud 
8  microphysics. 

Previous  studies  using  an  SE  model  (Giraldo  and  Restelli  2008)  have  used  a  fixed  set  of 
io  control  parameters  which  control  the  domain  decomposition:  the  number  of  elements  in  the 
horizontal  and  vertical  directions,  hx  and  hz,  respectively,  and  the  polynomial  order  p  (more 
12  details  in  section  2).  The  purpose  of  this  paper  is  to  assess  the  strengths  and  weaknesses  of  a 
SE  model  through  the  systematic  exploration  of  the  parameter  space  and  validate  simulation 
14  results  of  two  idealized  mesoscale  phenomena:  a  linear,  hydrostatic  mountain  wave  and  mid¬ 
latitude  squall  line. 

i6  In  the  first  part  of  this  paper,  we  use  an  analytic  solution  for  a  hydrostatic  mountain 
wave  to  validate  the  numerical  simulations  with  the  SE  model.  We  address  the  following 
is  questions:  1)  What  is  the  range  of  h  p  parameters  that  give  adequate  results?  2)  How 
computationally  expensive  is  the  SE  model  compared  to  a  typical  finite-difference  model? 
20  3)  How  quickly  does  the  SE  model  converge  to  the  final  solution? 

In  the  second  part,  we  assess  the  ability  of  the  SE  model  to  properly  simulate  the  squall 
22  line.  Since  there  is  no  analytic  solution,  the  typical  approach  with  numerical  simulations  is  to 
increase  both  spatial  and  temporal  resolution  until  convergence  is  achieved.  This  approach 
24  might  not  be  viable  for  atmospheric  convection  (e.g.  Weisman  et  al.  (1997),  Bryan  et  al. 

(2003)).  Therefore  we  focus  on  the  simulation  of  important  characteristics  (e.g.  cloud  forma- 
26  tion,  precipitation  initiation,  longevity  of  the  storm  system)  of  a  squall  line  and  integrated 
quantities  (e.g.  average  precipitation  accumulation,  maximum  vertical  velocity),  across  the 
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parameter  space. 


The  structure  of  this  paper  is  as  follows:  the  model  is  described  in  Section  2;  in  Section  3 
details  of  the  experiment  setup  are  given  followed  by  discussion  of  results  in  Section  4.  The 
conclusions  are  presented  in  Section  5. 


2.  Model  Description 

a.  Governing  Equations 

The  governing  equations  for  the  compressible,  nonhvdrostatic  numerical  model  of  the 
moist  atmosphere  are 

^  +  (Po  +  pOv'w  +  w^  +  w‘  Vp'  =  0  (1) 

du  1  o' 

—  +  u-Vu+ - —Vp'  +  gk( - — -0.61  q'v  +  qc  +  qr)  =  pV2-u  (2) 

ot  po  +  p'  po  +  p' 

80' 

-gf  +  n-Wv  =  SK+,j.V%,  (3) 

where  air  density  (p),  pressure  (p).  potential  temperature  (9)  water  vapor  mixing  ratio  (qv) 
are  separated  into  the  vertically  varying,  hydrostatically  balanced  base  state  (subscript  0) 
and  perturbation  (denoted  by  prime):  p  =  po(z)  +  p'(x,t).  The  wind  vector  has  a  horizontal 
and  vertical  component  u  =  (u(x,  t),  w(x,  t))T,  k  is  unit  vector  in  vertical  and  g=9.81  m 
s^2  is  the  acceleration  of  gravity. 

Moisture  related  variables  are  mixing  ratios  of  water  vapor  (qv),  cloud  water  (qc)  and 
rain  water  (qr).  They  are  predicted  according  to  a  simple  microphysical  mechanism  for  a 
warm  cloud,  that  does  not  include  ice,  snow  and  graupel,  as  follows 

d^  =  -Cc  +  Ec  +  Er  +  pV2qv 
at 

d^  =  Cc-Ec-Ac-Kc  +  pV2qc  (4) 

at 

——  =  Ac  +  Kc  —  Er  +  Fr  +  p\7~qr, 
at 

where  Cc  is  condensation  of  cloud  water,  Ec  is  evaporation  of  cloud  water,  Er  is  evaporation 


of  rain  water,  Ac  is  autoconversion  of  cloud  water  into  rain  water,  Kc  is  the  collection  of 
2  cloud  water  and  Fr  is  the  sedimentation  of  rain  drops  in  the  air  parcel  (Houze  1993).  For  a 
detailed  description  of  each  parameterized  process  see  Klemp  and  Wilhelmson  (1978). 

4  The  thermodynamic  equation  involves  a  source/sink  term  (Sgv)  that  describes  latent  heat 
release/uptake  during  phase  changes  of  moisture  variables.  The  momentum,  thermodynamic 
e  and  moisture  equations  involve  a  diffusion  term.  The  diffusivity  coefficient,  /x  =  200  m2s-1, 
represents  artificial  viscosity  terms  used  for  numerical  stability.  Note  that  the  diffusion  is  ap- 
s  plied  only  for  the  moist  sumulations.  The  pressure  ( p )  and  the  virtual  potential  temperature 
( 6V )  are  defined  as 

p  =  pRdT  (1  +  0. 61qv)  (5) 

R, i 

ev  =  (1  +  0.61  qv)T  (^j  ^  =  (1  +  0.61  qv)Q,  (6) 

io  with  T  being  the  air  temperature,  poo  =  105  Pa,  reference  air  pressure,  Rd  =  287  J  kg-1K_1| 
a  specific  gas  constant  of  dry  air  and  cp  =  1004  J  kg-1  specific  heat  at  constant  pressure  for 
12  dry  air. 

b.  Numerical  model 

14  Using  the  SE  model,  the  computational  domain  U  is  decomposed  into  Ne  nonoverlapping 
quadrilateral  elements  (Fig.  1,  left  panel) 

Ne 

u  =  |J  ne.  (7) 

e=l 

i6  Generally,  the  elements  do  not  need  to  be  quadrilateral  and  structured.  A  mapping  from  the 
global-domain  coordinate  system  x  =  (x,  z)  onto  the  element-local  (f2e)  coordinate  system 
is  £  =  (£,77)  is  described  by  an  element-specific  Jacobian  J  =  ||,  where  the  local  coordinates 
satisfy  (£,77)  G  [ — 1,  l]2  (Fig  1,  upper  right  panel). 

20  The  local  element-wise  solution  of  each  variable  /  can  be  discretized  using  ATh  order 
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polynomial  basis 


K 

fd(£,t)  =  J2MS)h(t),  /f  =  (Af  +  l)2,  (8) 

k= 1 

where  fd  is  a  discrete  representation,  ^  are  expansion  functions,  fk  are  expansion  coefficients  2 
and  N  +  1  is  the  number  of  expansion  functions  in  each  direction.  The  expansion  functions 
(Fig  1,  lower  right  panel)  are  constructed  as  4 

$k  =  hi(£(x))-hj(ri(x)),  i,j  =  l,...,N+l,  (9) 


where  ht  and  hj  are  Lagrange  polynomials 


MO 


1  (i-  0)^(0 

N(N  + 1)  (£-  0)^(0)’ 


(10) 


and  Pn( 0  are  Nth  order  Legendre  polynomials.  The  expansion  function  h%  is  zero  at  all  e 
nodal  points  except  0-  The  chosen  Legendre-Gauss-Lobato  (LGL)  grid  points  within  the 
elements  (0,M>  are  not  equally  spaced  (Fig  1,  upper  right  panel),  but  are  given  as  the  roots  s 
of 


(i  -  emo  =  0. 


(id 


The  LGL  points  with  the  associated  quadrature  weights  u ;(0) 

2 

/  /  1  \ 

MO)  = 


(12) 


N(N  +  1)  \PN(&), 

can  be  directly  used  for  the  Gaussian  quadrature,  approximating  integrals  over  a  local  ele¬ 
ment  Vte 

N+l 


u: 


/  f(x)dx  = 

The  governing  equations  to  be  solved  are  in  the  form 


f(€,v)J(£>v)d£dr}  ~  M0MM/(0>MIJ(0>MI-  (13) 

i,j=  1 


df_ 

dt 


+  F{f)  =  0. 


(14) 
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Substituting  for  the  discretized  solution  will  result  in  a  residual 

r)fd 

RUd)  =  -^  +  F(fi),  (15) 

2  which  can  be  minimized  by  various  methods.  In  the  Galerkin  method,  the  residual  is  or¬ 
thogonal  to  the  expansion  functions 

(R,fpk)  =  0,  k  =  1,...,(IV  +  1)2,  (16) 


where  the  Legendre  inner  product  (/,  g )  over  the  subdomain  Vte  is  defined  as 


{f,9)=  /  f{x)g{x)dx. 


(17) 


Combining  (8),  (15)  and  (16)  leads  to  a  system  of  differential  equations 


(Jv+i)2  , 

Et  aJk 

Ink  77"  - 


n=  1 


dt 


/(iV+l)2 

F  I  £  ^(0/n(*)  I  k  =  1,  .  .  .  ,  N  +  1. 


\  n= 1 


(18) 


The  orthogonality  of  expansion  functions  simplifies  the  calculation  of  the  mass  matrix  (I, 


nk  ) 


Ink  'Pm'tpk ) 


(19) 


The  right  hand  side  of  (18)  can  be  solved  using  the  Gaussian  quadrature.  The  spatial  deriva- 
s  fives  appearing  in  the  governing  equations  are  constructed  through  the  analytic  derivatives 
of  the  basis  functions  (for  brevity  only  one-dimensional  example  is  provided) 

8f  a/«s  9  =  (20) 


dx  d£  dx 


,fe= i 


\fc=i 


di 


dx 


io  c.  Time  integration 

The  left  hand  side  of  (18)  can  be  readily  integrated  in  time  with  a  desired  accuracy.  The 
12  nonuniform  spacing  of  the  nodal  points  can  impose  a  severe  constraint  on  a  time  step  when 
using  a  regular  explicit  time  integration  scheme.  For  example,  the  ratio  of  the  maximum  to 
i4  minimum  nodal  spacing  for  the  tenth  order  polynomial  expansion  functions  is  almost  five 
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and  the  maximum  time  step  required  for  numerical  stability  is  limited  by  the  minimum  nodal 
spacing.  The  terms  in  the  governing  equations  can  be  rewritten  in  a  compact  vector  form 


2 


dq 

Ot 


S(q), 


(21) 


where  q  =  (p',u1,9v,qv,qC}qr)T,  and  S  represents  all  terms  not  involving  time  derivatives. 
The  semi-implicit  time  integration  can  be  introduced  in  the  previous  equation  as 

=  {S(q)  -  \L(q)}  +  [AL(<j)l,  (22) 


4 


where  curly  and  square  braces  represent  explicit  and  implicit  integration,  respectively,  A  = 

{0, 1}  is  a  control  flag  to  invoke  the  implicit  integration  (A  =  1),  and  L  is  a  linear  approx-  e 

imation  of  S  that  contains  acoustic  and  gravity  waves.  The  moisture  related  variables  are 
not  responsible  for  either  fast  mode,  so  the  linearized  formulation  contains  only  the  first  four  a 
components  of  q 


L(q) 


_ 1  dp' 

Po  dx 


__L fV  _  nP 
po  dz  y 


PO 


dOyQ 


(23) 


Instead  of  solving  for  each  variable  separately,  they  are  combined  into  one  pseudo-Helmholtz  io 
equation  for  the  pressure  perturbation  (Schur  complement).  Upon  solving  for  the  pressure, 
each  of  the  prognostic  variables  can  be  solved  in  sequence  for  the  updated  variables  (Giraldo  12 
et  al.  2010).  The  time  integrator  is  the  second  order  backward  difference  method,  BDF2 
(Giraldo  2005).  It  is  used  in  a  semi-implicit  mode  permitting  longer  time  steps  compared  u 
to  fully  explicit  methods  with  equal  or  higher  order  of  accuracy  which  have  also  been  tested 
(e.g.  family  of  Runge-Ivutta  schemes)  16 

As  in  most  numerical  models  involving  moist  processes  solved  with  a  finite-difference 
scheme  in  the  vertical,  the  microphysics  computation  is  time-integrated  separately  to  allow  is 
a  time  step  adjustment  for  the  case  when  sedimentation  of  precipitable  water  is  too  fast. 
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Moist  processes  are  treated  in  a  column-wise  fashion,  descending  from  the  top  of  the  domain 
2  to  the  lowest  level,  moving  laterally  through  the  domain.  The  indexing  of  the  elements, 
and  all  the  loops  in  the  source  code  can  be  completely  unstructured  and  therefore  not 
4  readily  applicable  for  microphysics  calculations.  For  the  purposes  of  this  paper,  the  element¬ 
wise  thermodynamic  and  moist  variables  are  mapped  to  regular  two-dimensional  arrays 
e  suitable  for  column-wise  calculations.  Once  the  microphysics  computations  are  concluded, 
the  updated  variables  are  mapped  back  onto  their  local  elements.  As  such,  the  actual 
8  microphysical  processes  are  not  strictly  computed  within  the  semi-implicit  realm,  but  the 
advection  and  diffusion  of  the  moisture  related  variables  are. 

io  d.  Accuracy 

When  using  a  polynomial  expansion  basis,  one  frequently  refers  to  it  only  by  its  order.  It 
12  should  be  emphasized  that  this  is  not  the  same  ’order’  as  the  one  used  to  identify  the  leading 
term  of  the  error  when  using  finite-difference  schemes,  which  in  fact  describes  accuracy. 
14  Evaluation  of  Gaussian  quadrature  (RHS  of  (18))  over  N  +  1  quadrature  points,  will  be  exact 
to  machine  precision  as  long  as  the  polynomial  integrand  is  of  the  order  2  ■  (N  + 1)  —  3,  or  less 
i6  (Ivarniadakis  and  Sherman  2005).  Applying  the  SE  method  to  the  governing  equations  will 
result  in  an  inner  product  of  two  polynomials  of  the  same  (or  lesser)  degree.  For  example, 
is  the  product  u% subject  to  the  orthogonality  condition  is  of  the  3 N  —  1  order.  For  the 
exact  integration,  at  least  | N  +  1  points  are  needed,  while  only  N  +  1  are  available.  The 
20  integrand  is  subsampled  and  consequently  aliased.  To  eliminate  this  aliasing,  a  low-pass  filter 
is  applied,  but  not  directly  to  the  chosen  expansion  functions  because  they  are  nodal.  They 
22  are  transformed  into  modal  functions  first,  filtered  using  a  Boyd-Vandeven  filter  (Giraldo 
and  Rosmond  2004)  and  inversely  transformed  to  retrieve  a  filtered  set  of  nodal  expansion 
24  functions.  The  inexact  integration  has  a  very  minimal  impact  for  higher  order  polynomials 
(N  >  4)  (Giraldo  1998).  Note  that  an  exact  integration  could  be  achieved  by  using  a  separate 
26  set  of  quadrature  points,  but  the  accompanying  computational  cost  is  usually  prohibitive 
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due  to  the  mass  matrix  no  longer  being  diagonal.  The  errors  stemming  from  the  BDF2  time 
integration  are  of  second  order  accuracy. 


2 


3.  Setup  and  Initial  Conditions 

The  model  is  applied  in  a  two-dimensional  mode  with  a  horizontal  and  vertical  domain  of  4 
240  and  24  km,  respectively.  Due  to  the  irregular  spacing  of  nodal  points  within  an  element 
in  both  the  horizontal  and  vertical  directions,  a  nominal  resolution  is  introduced,  defined  as  e 
the  element’s  extent  divided  by  the  number  of  nodal  points  (minus  one)  in  that  direction. 
When  describing  a  simulation,  its  corresponding  nominal  horizontal  and  vertical  resolution  s 
are  provided.  The  discrepancy  between  the  actual  adjacent  nodal  point  spacing  and  the 
associated  nominal  resolution  increases  with  the  polynomial  order.  io 

The  desired  nominal  resolution  can  be  achieved  by  increasing  the  number  of  elements 
( h )  while  holding  the  polynomial  order  ( p )  constant  (’h-refinement’),  keeping  the  element  12 
number  constant  and  increasing  the  polynomial  order  (’p-refinement’),  or  varying  both.  The 
limiting  ’/?, -refinement’  case  is  a  finite-element  method  (high  /r,  low  p),  while  a  similar  limit  m 
for  the  ’p-refinement’  is  a  spectral  method  (one  element,  high  number  of  basis  functions). 

When  investigating  the  resolution  dependence  of  the  SE  model,  one  has  to  consider  ie 
exploration  of  the  phase  space  defined  by  both  parameters:  the  polynomial  order  and  number 
of  elements.  In  order  to  achieve  a  nominal  horizontal  resolution  representative  for  a  mesoscale  is 
model  we  choose  the  polynomial  order  to  vary  between  4  and  10  and  the  number  of  elements 
in  the  horizontal  direction  between  6  and  120.  The  resulting  nominal  grid  spacing  varied  20 
between  200  and  10000  m  in  91  simulations  overall  (see  Table  1  for  details). 

Note  that  the  same  refinement  applies  to  both  the  horizontal  and  vertical  directions.  It  22 
may  be  desirable  to  keep  the  nominal  vertical  resolution  constant  in  all  experiments,  such 
as  in  Weisman  et  al.  (1997),  to  focus  solely  on  the  effects  of  variations  in  the  horizontal  24 
resolution.  This  not  practical  to  two  reasons:  i)  the  polynomial  order  is  the  same  in  both 
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directions  in  the  current  version  of  the  model,  and  ii)  the  number  of  elements  can  be  an  integer 
2  number  only.  The  impact  of  varying  the  nominal  vertical  resolution  is  briefly  examined  in 
section  4  and  shown  to  have  significant  impact  as  long  as  the  nominal  vertical  resolution  is 
4  sufficient  to  adequately  resolve  the  squall  line  cold  pool.  The  ratio  of  nominal  horizontal  and 
vertical  resolution  is  thus  kept  in  the  same  range  (1:3-1:5)  in  all  simulations. 


e  a.  Dry  Experiments  -  Linear,  Hydrostatic  Mountain  Wave 


In  the  first  suite  of  experiments,  we  focus  on  the  case  of  a  linear  hydrostatic  gravity  wave 
generated  by  flow  over  small  topography.  The  topography  h(x)  is  defined  with  a  bell  shaped 
profile 


h  a2 

h(x)=  •“ 


(24) 


x2  +  a2  ’ 

where  hm  =  1  m  is  the  terrain  height  and  a  =  10.0  km  is  the  mountain  half-width.  In 
an  isothermal  atmosphere  (T0  =  250.0  K)  the  atmospheric  stability  is  constant  with  height 
(IV  =  -A —  =  0.0196  s-1).  By  choosing  an  appropriate  wind  speed  (u  =  20.0  m/s)  both 

y/CpTo 

conditions  for  hydrostacity  ( Na/u  S>  1)  and  linearity  (. Nhm/u  <C  1)  are  satisfied.  The 
computational  domain  is  240  km  wide  and  24  km  deep.  An  active  sponge  layer  at  the  lateral 
and  top  boundaries  helps  to  damp  reflections  from  the  domain  boundaries.  Since  the  width 
of  the  lateral  sponge  is  proportional  to  the  nominal  horizontal  grid  spacing,  the  domain 
extent  is  doubled  horizontally  for  all  cases  with  the  nominal  grid  spacing  greater  than  4  km. 
All  simulations  are  performed  using  the  semi-implicit  time  integrator.  Each  simulation  is 
integrated  for  12  hours  (nondimensional  time  ztt/a  =  86.4),  assuring  that  a  steady  state  is 
reached. 


b.  Moist  Experiments  -  Squall  Line 

22  The  initial  conditions  are  specified  by  a  synthetic  vertical  profile  (Fig.  2),  based  on  a 
typical  environment  for  midlatitude  squall  lines  and  used  in  several  previous  studies  (Rotunno 
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et  al.  1988).  It  features  increasing  moisture  in  the  lower  troposphere  and  a  fairly  moist  but 
unsaturated  air  mass  in  the  rest  of  the  troposphere.  The  air  is  weakly  stable  close  to  the  2 
ground  with  uniform  stability  (!V=0.01  s'1)  up  to  the  tropopause  at  12  km  where  the  stability 
increases  (IV=0.02  s'1).  A  low-level  wind  shear  is  added  to  promote  longevity  of  the  storm  4 
by  separating  the  storm  inflow  from  the  downdraft  created  by  precipitation.  In  addition,  if 
the  horizontal  component  of  vorticity  of  the  environmental  shear  is  approximately  balanced  s 
by  the  vorticity  of  the  opposite  sign,  created  by  the  density  current  of  the  outflow,  the  storm 
will  remain  quasi  stationary  (Rotunno  et  al.  1988).  The  topography  is  set  to  zero  for  the  8 
moist  experiments,  there  is  a  sponge  layer  at  the  top  of  the  domain,  identical  to  the  dry 
simulations,  but  at  the  lateral  boundaries  we  use  periodic  boundary  conditions.  The  main  10 
reason  for  choosing  periodic  lateral  boundary  conditions  is  to  evaluate  mass  conservation 
during  the  simulation  .  12 

The  triggering  mechanism  for  the  storm  evolution  is  a  warm  bubble  (Rotunno  et  al.  1988), 
centered  at  the  height  of  2  km,  inserted  at  the  initial  time.  The  temperature  perturbation  14 


is  defined  as 

/ 

9C  ■  cos2  r  <  rc 

A  9  =  < 

(25) 

0,  r  >  rc 

r  =  « 

\x-x 

(26) 

1 

1  xr  zr 

where  9C= 3.0  Iv,  xc=0,  zc=2,  xr=10, 

zr  1.5  and  rc=1.0  km. 

The  perturbation  reaches 

its  maximum  value  at  the  center  (xc,  zc)  and  decreases  radially  outward.  The  triggering 
mechanism  is  different  from  the  density  current  used  by  Weisman  et  al.  (1997).  Due  to  the  is 
periodic  boundary  conditions  used  in  our  simulations,  the  density  current  would  enter  the 
domain  from  the  upstream  and  cause  an  unwanted  secondary  line  of  storms.  20 

The  initial  positive  buoyancy  perturbation  initiates  air  parcel  ascent.  Once  they  reach 
the  level  of  free  convection,  the  lifting  continues  as  long  as  the  parcels  are  less  dense  than  the  22 
surroundings,  described  by  the  Convective  Available  Potential  Energy  (CAPE),  summarized 
in  Table  2.  Values  in  excess  of  2000  J  kg'1  suggest  the  possibility  of  a  strong  convective  24 
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activity. 


4.  Results 


a.  Dry  Experiments  -  Linear,  Hydrostatic  Mountain  Wave 


We  explore  the  h  —  p  parameter  space  through  analysis  of  the  inviscid  (i.e.  no  artificial 
viscosity  p,  —  0),  linear,  hydrostatic  mountain  wave  simulations,  for  which  an  analytic  so¬ 
lution  exists.  Instead  of  calculating  error  statistics  for  the  model  variables  separately,  the 
model  performance  is  assessed  by  calculating  a  second  order  quantity,  the  momentum  flux, 
as  a  function  of  height,  Mz,  and  compared  to  the  analytic  height-independent  solution  (Ma) 
(Durran  and  Ivlemp  1983) 


Mz 

Ma 


p(z)v!  (z)w(z)dx 


—  ^ p0Nuh 2m  =  —0.4285  kg  s  2, 


(27) 

(28) 


where  p0  =  1.3937  kg  m-3  is  the  air  density  at  the  surface  and  the  u  component  of  velocity 
is  decomposed  into  the  mean  state  and  perturbation  (u  —  u  +  u').  Due  to  the  small  terrain 
height  (hm),  the  quadratures  calculated  at  a  constant  height  z  (Mz)  or  at  a  constant  model 
level  k  (Mk)  are  essentially  the  same.  While  evaluating  the  integrals,  the  lateral  portions  of 
the  domain  with  an  active  sponge  are  omitted.  The  normalized  1 2  norm  is  calculated  using 
Mk  and  Ma  for  all  the  model  levels  from  the  ground  (k  =  1)  to  the  uppermost  level  not 
affected  by  the  top  sponge  (k  =  ks) 


h 


Eti(Ma)2 


2 


(29) 


Simulations  with  higher  nominal  resolution  (Ax  <  1  km)  have  the  smallest  error  statistics 
(right  portion  of  Fig.  3).  Cases  with  lower  nominal  horizontal  resolution  (Ax  >  5  km)  and 
polynomial  order  (p  <  6)  result  in  relatively  poor  l2  statistics  (lower  left  portion  of  Fig. 


3),  due  to  a  combined  effect  of  the  poorly  resolved  topography  and  error  introduced  with 
inexact  integration  (Eqn.  (13))  for  lower  polynomial  orders. 

A  smaller  subset  of  cases  (shaded  in  gray  in  Table  1)  is  further  analyzed  to  assess  the 
speed  of  convergence  and  computational  cost  as  a  function  of  the  polynomial  order  (p= 4, 
6,  8  and  10)  and  number  of  elements  (/?,),  reflected  in  the  nominal  horizontal  resolution 
(Aa;=0.5,  1.0,  2.0  and  3.0  km).  In  addition,  accuracy,  convergence  and  timing  comparisons 
for  a  finite-difference  (FD)  model  (fully  compressible,  nonhydrostatic  with  a  fourth  order 
horizontal  advection,  for  more  details  see  Durran  and  Ivlemp  (1983))  are  used  with  the  same 
domain  and  with  the  matching  horizontal  and  vertical  grid  spacings. 

A  reduction  of  the  l2  error  at  a  particular  point  on  Fig.  3  can  be  achieved  by  increasing  the 
polynomial  order  (p-refinement),  increasing  the  number  of  elements  (h  refinement)  or  both. 
Keeping  the  relatively  coarse  nominal  resolution  constant  (Ax=2.0,  3.0  km)  and  increasing 
p  yields  only  minor  error  reduction  (compare  dotted  and  dash-dotted  lines  of  various  shades 
of  gray  in  Fig.  4).  There  is  not  much  change  in  the  l2  error  past  t= 5  hours  regardless  of  the 
polynomial  order,  which  suggests  that  the  nominal  resolution  is  too  coarse  for  an  accurate 
solution  of  the  mountain  wave.  Starting  with  a  finer  resolution  (Ax=0.5  or  1.0  km)  the 
error  continues  to  decrease  and  it  is  at  least  an  order  of  magnitude  smaller  compared  to  the 
previously  analyzed  set  of  cases  with  coarser  resolution,  regardless  of  the  polynomial  order. 
Note  that  while  the  errors  for  the  Aa;=1.0  km  set  reach  minima  at  t=ll  h,  they  are  still 
decreasing  at  the  end  of  the  simulation  time  for  the  Aa;=0.5  km  set.  For  the  finite-difference 
model  (lines  with  diamonds  in  Fig.  4),  increasing  the  spatial  resolution  and  decreasing  the 
time  step  results  in  a  more  monotonic  error  reduction.  Even  at  the  highest  resolution  (solid 
line  with  diamonds  on  Fig.  4)  the  error,  dominated  by  the  lowest-order  truncation  term,  is 
still  at  least  an  order  of  magnitude  larger  compared  to  the  SE  model  (solid  lines  on  Fig.  4). 
Note  that  the  error  lines  for  the  simulations  with  p=10  indicate  the  error  is  still  decreasing 
after  12  hours  of  simulation,  previously  observed  by  Giraldo  and  Restelli  (2008). 

Improving  accuracy  by  using  more  elements  (better  resolution)  is  computationally  more 


costly.  Wallclock  time  for  the  SE  model  increases  approximately  by  an  order  of  magnitude,  as 
2  the  resolution  gets  refined  from  Ax=2.0  to  Ax=1.0  km  and  similarly  from  Ax=1.0  to  Ax=0.5 
km  (Fig.  5).  Variations  within  each  cluster  of  points  are  due  to  the  polynomial  order,  where 
4  the  lowest  order  is  the  least  expensive  (compare  matching  symbols  with  different  shades  of 
gray  in  Fig.  5).  The  finite-difference  model  is  computationally  less  expensive  compared  to 
e  the  SE  model,  when  comparing  timing  results  for  matching  nominal  spacing  (Ax)  for  the  SE 
model  with  the  constant  spacing  (Ax)  for  the  FD  model.  Note,  however,  that  if  we  change 
s  the  comparison  metric  to  a  desired  value  of  /2  error,  the  SE  model  is  faster.  Moreover,  at 
the  same  computational  cost,  the  /2  error  associated  with  the  SE  model  is  for  the  values  of 
io  Ax  <  1  km  at  least  an  order  of  magnitude  smaller  compared  to  the  FD  model.  The  error 
reduction  is  gradual  with  increasing  resolution  for  the  FD  model  (solid  line  with  triangles), 
12  while  the  major  error  reduction  for  the  SE  model  occurs  with  a  refinement  from  Ax=2.0  to 
Ax=1.0  km  (Fig.  5).  The  integration  time  steps  used  for  both  models  are  at  a  maximum 
14  allowed  from  a  numerical  stability  perspective. 

To  summarize  the  dry  experiments,  the  resolution  required  to  adequately  resolve  the 
i6  simulated  phenomenon  can  be  achieved  by  either  h  or  p  refinement.  At  a  fixed  nominal  reso¬ 
lution  the  error  is  almost  always  the  largest  for  the  lowest  polynomial  order,  p= 4,  represented 
is  by  black  lines  in  Fig.  4.  Our  recommendation  is  therefore  for  the  polynomial  order  to  be  at 
least  p= 6.  Higher  values  of  p  come  with  increasing  computational  cost,  perhaps  prohibitively 
20  expensive  for  p= 10,  with  the  best  ratio  of  accuracy  and  resources  spent  is  achieved  at  p= 8. 

The  number  of  operations  for  a  two-dimensional  SE  model  described  in  this  paper  is  on  the 
22  order  of  0(Ne  -p3),  with  Ne  being  the  total  number  of  elements  (a  product  of  number  of 
elements  in  the  horizontal  and  vertical).  As  the  resolution  refinement  scales  as  0(Ne  ■  p),  it 
24  is  computationally  more  feasible  to  increase  the  number  of  elements,  since  the  cost  increases 
cubically  with  p.  With  the  fixed  nominal  resolution,  the  ratio  of  the  most  expensive  (p=10) 
26  to  the  least  expensive  (p= 4)  simulation  is  2.5,  which  can  be  calculated  from  the  table  1  and 
confirmed  on  Fig.  5. 
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b.  Moist  Experiments  -  Squall  Line 


A  brief  synopsis  of  the  storm  evolution  is  based  on  a  simulation  with  a  typical  mesoscale  2 
resolution  with  Aa;=l  km  and  Ax:=0.2  km,  which  corresponds  to  a  case  with  p= 8  and 
h= 30.  After  the  initialization,  the  first  cloud  forms  at  around  t=900  s  (Fig.  6a),  the  cloud  4 
keeps  growing  with  warm  microphysical  processes  resulting  in  rain  formation,  which  starts 
accumulating  at  the  surface  at  around  t=1800  s  (Fig.  6b).  By  t=4800  s,  a  strong  cold  s 
pool  has  formed  at  near  the  rear  of  the  storm,  characterized  by  negative  equivalent  potential 
temperature  perturbation,  caused  by  evaporative  cooling  by  rain  water  and  downward  motion  8 

of  cooler  air  from  aloft  (Fig.  6c).  The  cold  pool  spreads  as  a  density  current  at  the  surface 
and  if  the  induced  shear  and  associated  horizontal  component  of  vorticitv  are  not  exactly  10 
balanced  by  the  ambient  shear,  the  squall  line  will  propagate.  New  updrafts  are  being  formed 
in  the  upshear  region  (ahead  of  the  location  of  the  original  initiation)  as  the  density  current  12 
initiates  forced  lifting  (secondary  triggering  mechanism),  consistent  with  a  broadening  region 
of  accumulated  precipitation  and  downshear  tilt  of  the  convective  tower  (Fig.  6d).  The  m 
subsequent  triggered  convection  is  generally  weaker  compared  to  the  initial  onset. 

We  start  examining  the  results  across  the  h-p  parameter  space  by  inspecting  simulations  ie 
with  the  same  nominal  resolution  Ax=l  km,  at  t=6000  s  (Figs  7a-f).  Overall,  the  cloud 
structure  (anvil  extent,  downshear  tilt  of  the  convective  tower),  the  cold  pool  intensity  and  is 
precipitation  amount  are  similar  among  the  simulations.  The  most  significant  difference 
among  the  simulations  is  the  spatial  distribution  of  the  rainfall  accumulation  and  related  20 
lateral  extent  of  the  cold  pool  beneath  the  cloud.  The  only  differences  in  the  setup  among 
the  cases  are  the  polynomial  order  p  and  the  number  of  elements  in  the  horizontal  direction  22 
h,  resulting  in  a  variable  nodal  spacing  where  the  ratio  of  the  maximum  to  minimum  nodal 
spacing  ranges  from  1.9  (p= 4,  Fig.  7a)  to  8.8  (p  20.  Fig.  7f).  Note  that  the  last  case  with  24 
high  value  of  p  is  not  described  in  Table  1.  Three  additional  experiments  are  designed  with 
p= 20.  The  number  of  horizontal  and  vertical  elements  is  6/3,  12/6  and  24/12,  resulting  in  26 
nominal  resolutions  of  2.0,  1.0  and  0.5  km,  respectively  (not  in  Table  1).  Despite  the  ratio  of 
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the  narrowest  to  the  widest  nodal  spacing  within  the  element  is  0(0.1)  (Table  3),  the  overall 
2  storm  is  still  well  resolved  (Fig.  7f).  The  rapidly  varying  nodal  spacing,  in  addition  to  large 
differences  between  the  widest  and  narrowest  nodal  distance,  does  not  result  in  preferred 
4  location  for  convection  or  “single-cell”  storms  or  updrafts. 

The  similarity  among  the  snapshots  of  the  squall  line  simulation  across  the  h-p  parame- 
e  ters  with  the  same  nominal  spatial  resolution  extends  the  robustness  of  the  SE  model  beyond 
the  dry,  dynamical  core  tests  shown  in  the  previous  section.  The  disagreement  in  the  total 
s  precipitation  accumulation  is  discussed  later  in  this  section. 

Adequately  resolved  storms  (cases  with  Ax  <  3  km)  undergo  similar  stages  of  devel- 
io  opment,  but  differ  in  the  accumulated  precipitation  amount,  as  shown  by  a  series  of  four 
simulations  with  the  same  polynomial  order  (p= 8).  The  nominal  resolution  starts  at  Ax  3 
12  km  and  is  progressively  reduced  by  factors  of  two  down  to  0.375  km  (Fig.  8).  The  simulation 
with  the  coarsest  nominal  resolution  has  an  excessive  amount  of  precipitation  with  an  over- 
14  all  cloud  outline  similar  to  the  shape  at  the  higher  resolution  (Fig.  8a).  As  the  resolution 
increases,  the  overall  precipitation  amount  decreases,  the  spatial  extent  of  the  cold  pool  is 
i6  reduced,  although  the  strength  is  comparable,  and  the  size  of  the  cloud  gets  smaller  (Figs. 8 
b-d). 

is  Without  an  existing  analytic  solution  for  comparison  purposes,  we  assess  the  moist  sim¬ 
ulation  using  metrics  appropriate  for  convective  events:  total  rain  accumulation,  maximum 
20  rain  rate  and  maximum  vertical  velocity.  Simulations  with  a  poorly  resolved  triggering 
thermal  bubble,  which  never  develop  any  convection  are  assigned  zeros  for  all  validation 
22  parameters.  These  simulations  are  clustered  in  the  lower  left  portion  of  the  h-p  parameter 
space  (Figs  9a,-c).  In  addition,  cases  filling  the  rest  of  the  void  region  share  in  common  the 
24  maximum  nodal  spacing  being  larger  than  4  km  (the  actual  limiting  contour  is  between  the 
Ax=2.o  and  3  km  contours).  This  latter  group  of  cases  grossly  overpredicts  the  precipitation 
26  (Fig  8a)  and  all  the  validation  parameters  are  assigned  zeros.  A  threshold  value  of  minimum 
grid  spacing  required  for  an  adequately  resolved  squall  line  is  similar  to  that  for  the  FD 


17 


models  (Weisman  et  al.  1997),  despite  the  difference  in  uniformity  of  grid  points  between 
the  SE  and  I'D  model. 

The  total  rain  accumulation  for  the  duration  of  the  simulation,  averaged  over  the  whole 
domain  (Fig.  9a)  indicates  a  decrease  with  increasing  h.  The  negative  correlation  between 
the  precipitation  accumulation  and  nominal  resolution  is  also  apparent  when  comparing  sim¬ 
ulations  with  the  same  p  (Figs.  8a-d).  Note  that  the  gradient  is  somewhat  independent  of 
the  polynomial  order.  The  reduction  of  the  rain  accumulation  is  consistent  with  findings 
of  Weisman  et  al.  (1997)  where  simulations  with  coarser  resolution  tended  to  exhibit  slower 
evolution,  stronger  storm  circulation  and  higher  overall  precipitation  amounts.  The  maxi¬ 
mum  rain  rate  (Fig.  9b)  is  reduced  with  increasing  h,  which  is  consistent  with  the  reduction 
of  the  total  precipitation  accumulation  at  higher  resolution.  A  comparison  of  the  maximum 
vertical  velocities  indicates  they  are  in  the  range  between  20  and  30  m  s_1  (Fig.  9c),  similar 
to  values  reported  by  Bryan  et  al.  (2006)  and  Weisman  and  Rotunno  (2004).  There  is  a  no¬ 
ticeable  trend  of  higher  maximum  vertical  velocities  (in  excess  of  30  m  s-1)  with  increasing 
nominal  resolution.  The  apparent  inconsistency  with  the  reversed  trend  in  vertical  velocities, 
compared  to  previously  observed  gradients  of  precipitation  accumulation  and  maximum  rain 
rate,  is  due  to  scaling  of  the  maximum  rain  rate  by  the  corresponding  nodal  spacing.  If 
similar  scaling  is  applied  to  the  vertical  velocities,  as  a  proxy  for  the  vertical  mass  flux,  there 
is  again  a  reduction  in  values  with  increasing  resolution  (not  shown). 

A  trend  that  can  be  recognized  from  Figs  9a,-c  suggests  that  results  are  more  dependent 
on  the  h  than  p  refinement  and  that  the  gradient  with  respect  to  h  is  consistent  for  all 
analyzed  quantities.  These  conclusions  differ  from  Weisman  et  al.  (1997),  but  in  their  study 
the  finest  resolution  is  1  km,  the  horizontal  grid  spacing  is  constant  and  more  importantly, 
the  sub-gridscale  mixing  is  parameterized. 

For  dry  simulations  of  a  density  current  with  increasing  resolution  (Straka  et  al.  1993),  the 
solutions  are  converging  towards  the  solution  obtained  with  the  finest  resolution.  Mixing 
has  a  strong  effect  on  the  overall  evolution  of  the  storm.  When  utilizing  a  sub-gridscale 


physical  mixing  that  scales  with  the  horizontal  resolution,  a  certain  degree  of  convergence  of 
2  solutions  can  be  expected  when  spanning  a  wide  range  of  horizontal  resolutions  (Weisman 
et  al.  1997).  If  the  resolution  is  progressively  refined,  the  results  might  become  different 
4  to  some  extent  as  documented  by  Bryan  et  al.  (2003)  and  in  this  paper.  We  hypothesize 
that  the  nonlinear  character  of  the  moist  processes  leads  to  this  behavior.  The  spatial  and 
e  temporal  distributions  of  buoyancy  perturbations  depend  on  localized  phase  changes,  which 
will  differ  among  simulations  with  different  nodal  point  distribution.  To  explore  this  issue 
8  further,  we  ran  an  additional  set  of  cases  (p= 8,  increasing  h)  based  on  the  setup  for  the  squall 
line  simulations,  except  with  no  moisture  at  the  initial  time.  We  calculated  power  spectra 
io  of  vertical  velocities,  averaged  in  time  and  height,  as  a  function  of  horizontal  wavelength. 

Since  the  model  data  is  on  a  nonuniformlv  spaced  grid,  it  is  resampled  with  the  horizontal 
12  spacing  that  approximates  the  narrowest  nodal  spacing.  If  the  above  hypothesis  holds,  the 
spectra  for  the  “dry”  squall  line  simulations  should  converge.  The  power  spectra  peaks  are 
14  at  the  same  wavelength  and  the  spectra  width  do  not  change  with  the  resolution  (Fig.  10b), 
except  for  the  case  with  Ax=3  km,  which  poorly  resolves  the  initial  thermal  bubble.  For  the 
i6  original  squall  line  simulations,  there  is  a  broadening  of  the  power  spectra  and  a  shifting  of 
the  maxima  towards  the  shorter  wavelengths  with  increasing  nominal  resolution  (Fig.  10a). 
is  Whether  this  trend  continues  or  if  the  spectra  collapse  with  further  resolution  refinement 
is  beyond  the  scope  of  this  research.  In  a  separate  subset  of  experiments  (p= 8,  increasing 
20  h)  with  no  latent  heat  release  or  uptake  permitted,  the  cloud  shapes  are  almost  exactly  the 
same  independent  of  the  spatial  resolution. 

22  The  average  time  between  sequential  discrete  updrafts  is  determined  by  local  maxima 
in  positive  vertical  velocities.  The  time  is  well  within  the  documented  range  of  Rotunno 
24  et  al.  (1988),  corresponding  to  their  “optimal  state”,  except  when  the  nominal  resolution  is 
less  than  1  km  the  average  time  becomes  longer,  because  the  subsequent  convective  cells 
26  take  longer  time  to  form.  In  addition,  the  time  between  the  initial  storm  triggering  and 
rain  reaching  the  ground  is  consistent  with  findings  in  the  literature  (Weisman  et  al.  1997) 
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throughout  the  parameter  space  (not  shown),  suggesting  that  the  choice  of  h  p  parameters 
does  not  affect  the  storm  triggering  by  the  initial  buoyancy  perturbation.  2 

One  of  the  concerns  when  simulating  severe  convection  using  a  variable  grid  is  develop¬ 
ment  of  preferential  locations  for  convection,  manifested  by  extrema  in  vertical  velocities.  4 
To  test  the  uniformity  of  the  spatial  distribution  of  vertical  momentum,  we  combined  the 
vertical  velocities  at  each  cell  into  0.5  m  s-1  bins  in  the  range  [4.75,20.75[  m  s-1,  for  all  s 
the  available  output  times.  Next,  bins  of  all  vertical  cells  at  a  fixed  horizontal  distance  are 
combined,  resulting  in  a  two-dimensional  histogram  revealing  a  spatial  distribution  of  oc-  a 
currence  for  a  particular  vertical  velocity  bin  (not  shown).  Most  of  the  motions  with  higher 
absolute  vertical  velocities  are  occurring  in  the  eastern  part  of  the  domain,  as  expected,  10 
where  the  squall  line  slowly  propagates.  There  is  no  visible  evidence  of  convection  triggering 
at  preferred  locations  (narrowest  nodal  spacing  next  to  the  element  boundaries).  Further-  12 
more,  bins  of  all  the  cells  with  the  same  horizontal  dimension  are  combined  and  normalized 
to  obtain  a  relative  frequency  histogram  as  a  function  of  the  vertical  velocity  (not  shown).  14 
If  the  convection  is  indeed  taking  place  closer  to  the  element  boundaries  where  the  nodal 
spacing  is  at  a  minimum,  this  would  manifest  itself  in  the  histogram  by  a  higher  (lower)  ie 
relative  frequency  of  the  narrower  (wider)  cells,  which  is  not  the  case. 

As  mentioned  in  section  3,  the  number  of  elements  can  be  independently  set  in  both  ia 
directions,  changing  the  respective  nominal  resolution.  We  designed  a  small  subset  of  four 
experiments  based  on  a  case  with  p=10  and  Aa;=l  km  to  assess  the  effect  of  varying  nominal  20 
vertical  resolution  with  the  nominal  horizontal  grid  spacing  held  constant.  The  number  of 
elements  in  the  vertical  direction  is  4,  10,  20  and  40,  resulting  in  nominal  vertical  resolution  of  22 
600,  240,  120  and  60  m,  respectively.  The  horizontal  location  of  the  most  intense  convection, 
the  magnitude  of  the  maximum  updraft  and  the  overall  precipitation  accumulation  are  not  24 
sensitive  to  the  vertical  resolution. 

In  a  series  of  additional  tests,  sensitivity  to  domain  length,  symmetry,  wind  shear  and  26 
viscosity  are  investigated.  The  choice  of  periodic  boundary  conditions  does  not  have  a  sig- 
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nificant  impact  on  the  results,  when  compared  to  simulations  with  triple  the  domain  length. 
2  With  no  ambiental  wind  shear,  a  symmetric  storm  cloud  is  expected,  but  an  asymmetry  can 
develop  if  the  initial  thermal  bubble  perturbation  is  not  centered  exactly  over  the  symmetric 
4  nodal  points.  If  the  wind  shear  is  too  strong,  no  storm  develops,  similar  to  findings  of  Ro- 
tunno  et  al.  (1988)  and  Weisman  et  al.  (1988).  Small  values  of  viscosity  (p<5  m2  s-1)  can 
e  lead  to  numerical  instabilities,  while  large  values  (//> 750  m2  s_1)  inhibit  convective  activity. 


5.  Conclusions 

a  In  this  paper  we  examine  the  characteristics  of  a  two-dimensional  spectral  element  (SE) 
model  for  dry  and  moist  mesoscale  atmospheric  test  cases:  a  linear,  hydrostatic  mountain 
io  wave  and  a  squall  line,  respectively. 

There  are  two  parameters  that  control  the  setup  of  the  SE  model:  the  number  of  ele- 
12  ments  into  which  the  computational  domain  is  subdivided,  and  polynomial  order  of  the  basis 
functions  ( p ),  which  determines  the  number  of  nodal  points  within  the  element.  The  spatial 
14  resolution  for  the  SE  model  is  determined  by  the  choice  of  the  two  parameters  with  ranges 
from  4  to  10  ( p )  and  6  to  120  (number  of  elements  in  horizontal,  /?,),  resulting  in  the  average 
i6  horizontal  (vertical)  resolution  ranging  from  200  (40)  to  10000  (1500)  m,  and  a  total  of  91 
simulations  spanning  the  h-p  parameter  space. 

is  For  the  linear  hydrostatic  mountain  wave  case, an  analytic  solutionis  used  to  validate  the 
model  performance.  Generally,  cases  with  the  nominal  resolution  less  than  2  km  yield  the 
20  best  results,  with  no  significant  gain  in  accuracy  if  the  resolution  is  refined  beyond  1  km.  The 
least  skillful  results  are  attributed  to  coarse  resolution,  not  sufficient  to  resolve  the  mountain 
22  barrier,  and  to  the  low  polynomial  order,  which  contributes  to  the  error  when  using  the 
inexact  integration.  Simulations  with  coarser  nominal  resolution  converge  faster  towards  the 
24  steady  state  solution,  but  with  larger  error.  In  addition,  the  SE  model  results  are  compared 
to  solutions  obtained  by  a  finite-difference  (FD)  model  with  matching  spatial  resolutions, 
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for  accuracy  and  timing  purposes.  The  error  for  the  FD  simulations  monotonically  decreases 
with  refined  spacing,  but  even  with  the  finest  grid  spacing  (0.5  km),  the  error  is  an  order  of  2 
magnitude  larger  compared  to  the  fine  resolution  cluster  of  the  SE  model.  At  a  given  resolu¬ 
tion,  matching  the  nominal  spacing  (Ax)  of  the  SE  model  with  the  constant  spacing  (Ax)  of  4 
the  FD  model,  the  SE  model  is  approximately  an  order  of  magnitude  more  computationally 
expensive  than  the  FD  model.  The  situation  is  reversed,  if  a  specified  error  is  desired  -  the  s 
SE  models  is  less  expensive.  Moreover,  the  error  of  the  SE  model  for  the  nominal  spacing 
Ax  <  1  km  is  an  order  of  magnitude  lower  compared  to  the  FD  model  at  the  same  reso-  a 
lution.  The  computational  cost  as  a  function  of  the  grid  spacing  and  associated  time  step 
increases  almost  uniformly  for  the  FD  model,  while  there  is  almost  no  improvement  in  error  10 
with  associated  computational  cost  when  increasing  the  nominal  resolution  from  3  to  2  km 
or  from  1  to  0.5  km  for  the  SE  model.  The  best  improvement  occurs  when  the  resolution  is  12 
refined  from  2  to  1  km. 

Simulations  that  adequately  resolve  the  initial  warm  bubble  perturbation  for  the  squall  14 
line  case,  successfully  simulate  the  upscale  transition  from  a  local,  isolated  convective  cell 
into  a  mesoscale,  organized  storm  system.  The  results  of  the  main  set  of  moist  experiments  ie 
and  additional  sensitivity  tests  suggest  the  overall  ability  of  the  SE  model  to  adequately 
simulate  the  squall  line.  Increasing  the  nominal  resolution  below  1  km  leads  to  some  dif-  is 
ferences.  Qualitatively,  the  cloud  shapes  are  very  similar,  but  simulations  with  the  finest 
nominal  resolution  tend  to  produce  stronger  maximum  vertical  velocities  with  more  local-  20 
ized  and  reduced  precipitation  accumulation.  We  hypothesize  and  offer  evidence  that  this 
behavior  can  be  explained  by  a  nonlinear  nature  of  latent  heating  and  localization  of  buoy-  22 
ancy  sources.  A  comparison  of  averaged  power  spectra  of  vertical  velocity  for  the  original 
squall  line  simulations  and  a  modified  set  with  no  initial  moisture  indicates  shifting  of  the  24 
power  spectra  toward  smaller  scales  for  the  moist  cases,  when  the  resolution  is  refined.  How 
would  a  continuing  resolution  refinement  affect  power  spectra  for  the  original  squall  line  26 
remains  an  open  question.  At  this  point  simulations  with  very  high  spatial  resolution  are 
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computationally  too  expensive  to  run  with  a  serial  code. 

2  The  SE  model  supports  both  structured  and  unstructured  grids.  The  accuracy  can  be 
adjusted  with  the  same  code  by  choosing  the  control  parameters  ( h  and  p).  Unlike  the  nu- 
4  merical  models  that  use  the  terrain-following  vertical  coordinate,  the  SE  model  can  handle 
complex  topographical  features  with  extreme  slope  angles,  such  as  in  urban  environments, 
e  For  all  these  advantages,  there  is  a  price  to  pay.  The  source  code  is  generally  less  straight¬ 
forward  to  understand  compared  to  the  source  code  of  the  FD  models.  As  mentioned  earlier, 
s  the  SE  models  are  computationally  more  expensive  compared  to  the  FD  models,  when  used 
at  the  same  spatial  resolution. 

io  A  recommended  subspace  of  the  h  p  parameter  space  depends  on  a  compromise  among 
acceptable  error,  computational  cost,  and  required  resolution  to  resolve  the  feature  of  choice. 
12  Based  on  our  results  for  inviscid,  dry,  and  viscous,  moist  simulations  of  the  mesoscale 
phenomena,  which  are  dimensionally  similar,  the  nominal  resolution  should  be  within  the 
14  Aa:=0.5-2  km  range  and  the  polynomial  order  in  the  range  p=5-10.  This  study  is  to  our 
knowledge  the  first  attempt  to  systematically  map  the  h  p  parameter  space  for  using  the  SE 
i6  model  in  mesoscale  atmospheric  modeling. 

The  results  are  certainly  encouraging  enough  to  warrant  further  investigations  in  using 
is  the  SE  model  for  more  realistic  mesoscale  atmospheric  modeling  scenarios.  The  model  in 
its  dry  and  inviscid  form  is  currently  being  tested  in  three  dimensions  and  on  massively- 
20  parallel  computers.  This  will  allow  us  to  extend  the  parameter  space  to  include  very  fine 
spatial  resolutions  which  are  prohibitively  expensive  in  the  serial  mode.  In  the  future,  we 
22  plan  to  adapt  the  microphysics  scheme  to  three  dimensions,  expand  it  to  include  the  ice 
phase  and  implement  a  sub-gridscale  mixing  parameterization. 
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Table 


Height  (m) 

0 

500 

1000 

1500 

2000 

CAPE  (J/kg) 

2383 

2426 

2781 

1968 

1133 

Initial  air-parcel  heights  and  corresponding  CAPE  values. 
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p 

^%max  /  ^%min 

^%max  / 

^%min/ 

4 

1.90 

1.31 

0.69 

5 

2.43 

1.43 

0.59 

6 

2.76 

1.41 

0.51 

7 

3.26 

1.47 

0.45 

8 

3.62 

1.45 

0.40 

9 

4.11 

1.49 

0.36 

10 

4.48 

1.48 

0.33 

20 

8.77 

1.53 

0.17 

Table  3.  Polynomial  orders  (p)  and  associated  nodal  spacing  ratios. 
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a)  Total  rain  accumulation  (in  mm)  in  6  hours  averaged  over  the  whole  do¬ 
main,  b)  Maximum  rain  rate  (in  kg  m-1  s_1)  and  c)  Maximum  vertical  ve¬ 
locities  (in  m  s~x),  as  a  function  of  the  polynomial  order  (p)  and  number  of 
elements  in  horizontal  (h)  for  the  squall  line  simulations.  Curved  black  lines 
represent  constant  nominal  horizontal  resolution  (Ax,  or  constant  number  of 
nodal  points  in  the  horizontal  direction).  44 

10  Power  spectra  for  simulations  with  p= 8  and  varying  nominal  resolutions:  Ax=3.0 
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dark  grey  line)  and  Ax=0.375  km  (thinnest  black  line).  Panel  a)  is  for  the 
control  squall  line  simulations  and  panel  b)  is  for  the  “dry”  squall  line  (see 
text  for  further  explanation).  The  spectra  are  averaged  over  height  (0-12  km, 
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Fig.  1.  An  example  of  element  decomposition  of  the  computational  domain  for  a  case  with 
polynomial  order,  p= 6,  number  of  elements  in  horizontal  h= 20  and  10  in  vertical  (left  panel). 
Distribution  of  LGL  nodal  points  within  a  canonical  element  with  p= 6  (top  right  panel). 
Basis  functions  -0i , . . . ,  -0?  for  p= 6  (bottom  right  panel). 


36 


16 


15 

14 


13 

12 


11 


bjO 

•  f-H 

<d 

ffi 


8 

7 

6 


5 


4 

3 

2 

1 

0 


-12  0 

Temperature  (°C)  Wind  Speed  (m/s) 


Fig.  2.  The  synthetic  sounding  used  to  initialize  the  model.  Temperature  and  dew  point 
temperature  are  represented  by  a  thicker,  solid  black  and  dashed,  grey  line,  respectively. 
The  wind  speed  profile  is  given  on  the  right  panel. 
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^-refinement  (polynomial  order) 


Fig.  4.  Time  evolution  of  the  normalized  1 2  norm  of  momentum  flux  for  the  two-dimensional 
linear  hydrostatic  mountain  wave  simulations  (output  every  hour,  starting  at  t  =  1  h). 
Results  based  on  simulations  with  the  same  horizontal  resolution  are  grouped  by  a  line  type: 
dotted  (3.0  km), dash-dotted  (2.0  km), dashed  (1.0  km),  and  solid  (0.5  km). Shades  of  gray 
depict  the  polynomial  order  of  basis  functions:  black  ( p  =  4),  dark  gray  ( p  =  6), medium 
gray  ( p  =  8), and  light  gray  ( p  =  10). In  addition,  results  obtained  with  the  finite-difference 
model  (lightest  gray)  have  added  diamonds,  which  also  replace  dots  in  corresponding  line 
styles. 


39 


Normalized  time 


Fig.  5.  Normalized  /2  norm  of  momentum  flux  for  the  two-dimensional  linear  hydrostatic 
mountain  wave  simulations  as  a  function  of  normalized  computational  time.  Results  obtained 
with  the  SE  model  are:  solid  black  line  (p=4),  dashed  dark  grey  (p= 6),  dot-dashed  medium 
gray  (p=S)  and  dotted  lightgrey  (p=10).  The  lightest  gray  line  with  triangles  is  for  the 
finite-difference  model.  Simulations  with  the  same  nominal  resolutions  are  grouped  together 
and  represented  with  small  circles  (0.5  km),  circles  with  wide  rings  (1.0  km),  circles  with 
thick  inner  and  thin  outer  rings  (2.0  km)  and  circles  with  two  thick  rings  (3.0  km). 
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Fig.  6.  Vertical  cross  sections  of  the  squall  line  evolution  as  depicted  by  a  simulation  with 
p= 8,  h= 30,  Ax  1  km  Az=0.2  km  at  a)  900,  b)  1800,  c)  4800  and  d)  9000  s.  Filled  contours 
represent  equivalent  potential  temperature  perturbation  (c.i.  3  K),  positive  values  with  dark 
and  negative  values  with  light  contour  lines.  The  interval  centered  around  zero  ([—3,3])  is 
omitted).  The  cloud  water  mixing  ratio  ( qc  =  10~5)  in  thick  black  line  represents  the  outline 
of  the  cloud.  The  bottom  portion  of  each  panel  shows  rain  water  accumulation  as  a  function 
of  distance.  Only  a  smaller  subset  of  the  full  domain  is  shown  to  emphasize  the  details. 
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Fig.  7.  Same  as  Fig.  6,  but  at  time  £=6000  s  and  for  cases  a)  p  4,  h= 60,  b)  p= 5,  h= 48, 
c)  p= 6,  h= 40,  d)  p= 8,  h= 30,  e)  p=10,  /?.=24  and  f)  p=20,  h=12.  All  cases  have  the  same 
nominal  horizontal  resolution  Ax=l  km  and  time  step  A£=0.25  s. 
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Fig.  10.  Power  spectra  for  simulations  with  p= 8  and  varying  nominal  resolutions:Ax=3.0 
km  (thickest  light  grey  line),  Ax=1.5  km  (thick  grey  line),  Ax=0.7-5  km  (thin  dark  grey  line) 
and  Ax=0.375  km  (thinnest  black  line).  Panel  a)  is  for  the  control  squall  line  simulations 
and  panel  b)  is  for  the  “dry”  squall  line  (see  text  for  further  explanation).  The  spectra 
are  averaged  over  height  (0-12  km,  with  0.5  km  increment)  and  time  (0-4  h,  with  300  s 
increment). 
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